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We have conducted A-body simulations of the growth of Milky Way-sized halos in cold and 
warm dark matter cosmologies. The number of dark matter satellites in our simulated Milky Ways 
decreases with decreasing mass of the dark matter particle. Assuming that the number of dark 
matter satellites exceeds or equals the number of observed satellites of the Milky Way we derive 
lower limits on the dark matter particle mass. We find with 95% confidence m s > 13.3 keV for 
a sterile neutrino produced by the Dodelson and Widrow mechanism, m s > 8.9 keV for the Shi 
and Fuller mechanism, m s > 3.0 keV for the Higgs decay mechanism, and ttiwdm > 2.3 keV for a 
thermal dark matter particle. The recent discovery of many new dark matter dominated satellites 
of the Milky Way in the Sloan Digital Sky Survey allows us to set lower limits comparable to 
constraints from the complementary methods of Lyman-a forest modeling and X-ray observations 
of the unresolved cosmic X-ray background and of dark matter halos from dwarf galaxy to cluster 
scales. Future surveys like LSST, DES, PanSTARRS, and SkyMapper have the potential to discover 
many more satellites and further improve constraints on the dark matter particle mass. 



I. INTRODUCTION 

Cold dark matter (CDM) is extremely successful at de- 
scribing the large scale features of matter distribution in 
the Universe but has problems on small scales. Below 
the Mpc scale CDM predicts numbers of satellite galax- 
ies for Milky Way-sized halos about an order of mag- 
nitude in excess of the number observed. This is the 
'missing satellites' problem 0, One proposed solution 
is that, due to feedback mechanisms, some dark matter 
satellites do not form stars and are nonluminous dark 
halos [3- 7]. Another solution is the power spectrum of 
density fluctuations may be truncated which may arise 
if the dark matter is 'warm' (particle mass ~ 1 keV) 
instead of 'cold' (particle mass ~ 1 GeV). Warm dark 
matter (WDM) particles decouple from the other parti- 
cle species in the early Universe with relativistic velocities 
and only become nonrelativistic when about a Galactic 
mass (~ 10 12 M Q ) is within the horizon. Streaming mo- 
tions while the particles are still relativistic can erase 
density fluctuations on sub-Galactic scales and reduce 
the number of satellites. WDM models have been stud- 
ied by a number of authors m relation to the miss- 
ing satellites problem and other issues with CDM such 
as the apparent density cores in spiral and dwarf galaxies 

A^-body simulations of WDM cosmologies are compli- 
cated by numerical artifacts produced by the discrete 
sampling of the gravitational potential with a finite num- 
ber of particles (see [HJ for a review). Matter perturba- 
tions collapse and form filaments with nonphysical ha- 
los separated by a distance equal to the mean particle 
spacing (see Fig. [1} [U,[2^]. These halos are numerical 
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artifacts. The ability of these halos to survive disruption 
as they accrete from filaments onto Milky Way-sized ha- 
los has not been studied and they may contaminate the 
satellite abundances and distributions in WDM simula- 
tions. 



FIG. 1. Nonphysical halos formed along a filament and ac- 
creting onto a larger halo at z — 1 in a WDM simulation 
{raw dm = 1 keV). These halos are numerical artifacts. 

In the past few years, 16 new dwarf spheroidal galax- 
ies have been discovered in the Sloan Digital Sky Survey 
(SDSS) HH (see Table 3 and references therein). Af- 
ter correcting for completeness the estimated number of 
Milky Way (MW) satellites is > 60 (see Sec. lmC|) . These 
new dwarfs have low luminosities, low surface bright- 
nesses, and most appear to be dark matter dominated. 
Since the number of dark matter halos must be greater 
than or equal to the number of observed satellites, the 
new data from the SDSS may provide improved limits 
on the mass of the dark matter particle independent of 
complementary techniques. 

Motivated by the recent increase in the number of ob- 
served Milky Way satellites, we have performed new sim- 
ulations of the growth of Milky Way-like galaxies in CDM 
and WDM cosmologies for a variety of WDM particle 
masses. Our goal is to constrain the dark matter parti- 
cle mass by comparing the number of satellite halos in 
the simulated Milky Ways to the observed number of lu- 
minous satellites for the actual Milky Way. Maccio and 
Fontanot fl4j combined ./V-body simulations with semi- 
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analytic models of galaxy formation to compare the sim- 
ulated and observed Milky Way satellite luminosity func- 
tions for CDM and WDM cosmologies. In this work, we 
do not make any assumptions on how we populate dark 
matter halos with luminous galaxies. We simply impose 
that the number of observed satellites is less than or equal 
to the number of dark matter halos for a range of Galac- 
tocentric radii. This guarantees a robust lower limit on 
the dark matter particle mass. 



II. SIMULATIONS 

All our simulations were conducted with the iV-body 
cosmological simulation code GADGET-2 [27j assuming 
dark matter only. We adopted values for cosmological pa- 
rameters from the third year release of the WMAP mis- 
sion |H, i^m, Ha, h, a 8 , n s ) = (0.238, 0.762, 0.73, 0.751, 
0.951) to facilitate comparison with the Via Lactea II 
simulation (29j . For each simulation set we produced a 
single realization of the density field in the same periodic, 
comoving volume but varied the power spectrum of fluc- 
tuations appropriate for CDM and WDM cosmologies. 
Our initial conditions were generated on a cubic lattice 
using the GRAFIC2 software package [30j| . The power 
spectra for CDM and WDM are given by 
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respectively, with the normalization of Pcdm determined 
by erg. For our first two sets of simulations (set A and B 
described below) we used the transfer function for CDM 
adiabatic fluctuations given by Bardeen et al. (BBKS) 
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where q = k/(Q m h 2 ). A potential problem with the 
BBKS transfer function is that it underestimates power 
on large scales. In the Appendix we investigate the ef- 
fect that our choice for the CDM transfer function may 
have on the number of Milky Way satellites. We run one 
of our simulations adop ting the transfer functions from 
Eisenstein and Hu |32| and we find that this does not 
affect our results on the number of satellites. We also 
ran additional CDM simulations (set C) using the trans- 
fer function calculated by the LINGER program in the 
GRAFIC2 package (f2t = 0.04 was used for calculating 
the effects of baryons on the matter transfer function). 
LINGER integrates the linearized equations of general 
relativity, the Boltzmann equation, and the fluid equa- 
tions governing the evolution of scalar metric perturba- 
tions, photons, neutrinos, baryons, and CDM. The mass 
and circular velocity functions for satellites are consistent 
across both transfer functions. 



Assuming the WDM to be a thermal particle, a particle 
that was in thermal equilibrium with the other particle 
species at the time of its decoupling, we used the trans- 
fer function valid for thermal particles given by Bode, 
Ostriker, and Turok [Toj]: 
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The parameter gx is the number of degrees of freedom 
for the WDM particle, conventionally set to the value 
for a light neutrino species: gx = 1.5. The parameter k 
is the spatial wavenumber in Mpc -1 and mwDM is the 
mass of the WDM particle in keV. 

A candidate for a thermal WDM particle is the grav- 
itino, the superpartner of the graviton in supcrsymmetry 
theories. The lightest stable particle (LSP) in supersym- 
metry theories is a natural dark matter candidate. If the 
scale where supersymmetry is spontaneously broken is 
< 10 6 GeV, as predicted by theories where supersymme- 
try breaking is mediated by gauge interactions, then the 
gravitino is likely to be the lightest stable particle and 
can have a mass reaching into the keV regime [33T ]. 

In general the dark matter particle may not have been 
in thermal equilibrium when it decoupled. This is the 
case for a sterile neutrino (see [H| and references therein) , 
a theoretical particle added to standard electroweak the- 
ory, the only matter it interacts with (except through 
gravity) is left-handed neutrinos. Sterile neutrinos have 
been proposed [35T - [42| as an explanation for the anoma- 
lous excess of oscillations observed between muon and 
electron neutrinos and antineutrinos (43Tj49| . There are 
several mechanisms by which sterile neutrinos can be pro- 
duced. In the standard mechanism proposed by Dodelson 
and Widrow (DW) [50|, sterile neutrinos are produced 
when oscillations convert some of the more familiar active 
neutrinos into the sterile variety. The amount produced 
depends on the sterile neutrino mass and the mixing an- 
gle but we will not consider such details here and when 
considering sterile neutrinos we simply assume they com- 
pose the entirety of the dark matter. The transfer func- 
tion for DW sterile neutrinos with mass m s is given by 

M 
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where v = 2.25, /i = 3.08, and 
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a = 0.1959 
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Viel et al. [52| give a scaling relationship between the 
mass of a thermal particle and the mass of the DW ster- 
ile neutrino for which the transfer functions are nearly 
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identical: 
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Other sterile neutrino production mechanisms include 
that of Shi and Fuller (SF) jH who showed the DW 
mechanism is altered in the presence of a universal lep- 
ton asymmetry where production can be enhanced by 
resonance effects. Sterile neutrinos can also be produced 
from decays of gauge-singlet Higgs bosons at the elec- 
troweak scale [54[. The momentum distribution of the 
sterile neutrinos depends on the production mechanism. 
In the absence of transfer function calculations we use 
the expressions in (3~H for the free streaming length and 
average momentum to derive approximate scaling fac- 
tors for the SF and Higgs produced sterile neutrinos: 
m D w/msF = 1-5, m DW /m Higgs = 4.5. 

There are also ways other than WDM to reduce small 
scale power. Broken-scale invariance inflation models [55T | 
have a cutoff length below which power is suppressed. 
Particle theories where the LSP dark matter particle 
arises from the decay of the next lightest supersymmet- 
ric particle (NLSP) can also suppress small scale power if 
the NLSP is charged and coupled to the photon-baryon 
plasma [56| or if the NLSP decay imparts a large velocity 
to the LSP [57J. Further possibilities include composite 



dark matter models where stable charged heavy leptons 
and quarks bind to helium nuclei by Coulomb attraction 
and can play the role of dark matter with suppression 
of small scale density fluctuations [58U641 ]. Our method 
can potentially be applied to constrain these models as 
well. However, in the rest of this paper we will not dis- 
cuss further the consequences that our work has on these 
theories. 

In our simulations we assume the dark matter is ther- 
mal and scale the results to the standard sterile neutrino 
mass using Eq. ([5]). Our initial conditions include par- 
ticle velocities due to the gravitational potential using 
the Zeldovich approximation but we do not add random 
thermal velocities appropriate for WDM to the simula- 
tion particles. Bode, Ostriker, and Turok argue that 
for warm particle masses greater than 1 keV thermal mo- 
tions are unimportant for halos on scales of a kiloparsec 
and above. Regardless, we expect thermal motions, if 
anything, would reduce the number of small mass halos 
and by not including thermal motions the mass limits 
derived from our simulations will be more conservative. 

We ran simulations for CDM and WDM cosmologies 
with particle masses of my/DM = 1, 2, 3, 4, and 5 keV 
(m s = 4.4, 11.0, 18.9, 27.8, 37.4 keV). Figure [2] shows 
the power spectra for these cosmologies along with the 
spectrum for an 11 keV standard sterile neutrino using 
Eq. (|HJ|. We ran two separate sets of simulations both 
consisting of a comoving cubic box 90 Mpc on a side. 
Set A consisted of 204 3 particles giving a 'coarse' par- 
ticle mass of 3.0 x 10 9 M Q and a force resolution of 
8.8 kpc (all our force resolutions were fixed in comov- 
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FIG. 2. Power spectra for our simulations. The dotted line is 
the power spectrum for an 11 keV standard sterile neutrino 
from Abazajian [5J. The neutrino spectrum is approximately 
the same as a 2 keV thermal particle, validating the scaling 
relation of Viel et al. [5^] . The vertical dashed lines indicate 
the lattice cell size in our high and low resolution refinement 
levels. 



ing coordinates) . We ran the HOP halo finding software 
65| at z = and identified Milky Way-sized halos with 
masses 1 — 2 x 10 12 A/ Q . Halos were examined visually, 
one was chosen that was at least several Mpc away from 
clusters and other large structures so as to be relatively 
isolated. Its particles were identified in the initial con- 
ditions and a cubic refinement level, 6.2 Mpc on a side, 
was placed on the region. For the refinement region in 
our low resolution simulations we used 11, 239, 424 (224 3 ) 
particles with mass and force resolutions of 7.3 x 10 5 M© 
and 550 pc, respectively. For CDM and WDM parti- 
cle masses of 1, 2, and 4 keV we also conducted higher 
resolution simulations with 89,915,392 (448 3 ) particles 
in the refinement region and mass and force resolutions 
of 9.2 x 10 4 M Q and 275 pc, respectively. The simulated 
Milky Way halo had a neighbor halo with mass 0.23Mmw 
at a distance of 700 kpc in the low resolution simulations. 
The real Milky Way has a massive neighbor in M31, the 
Andromeda galaxy (Mm31 ~ 1 — 3Mmw), at a distance 
~ 700 kpc. Being nonlinear and chaotic systems, small 
perturbations to the trajectories of dark matter halos can 
be amplified exponentially and in the higher resolution 
simulation this satellite is merging with the Milky Way 
at z = 0. Such a merger may disrupt the equilibrium 
of the halo and make it nonrepresentative of the actual 
Milky Way. The difference between the high and the 
low resolution simulations is significant and complicates 
the comparison between the resolutions; however, this 
merger is not a violation of the selection method used 
for the set C halos described below and we find excellent 
agreement across all simulation sets. 

The need to explore the scatter between the subhalo 
distributions of different realizations of Milky Way-type 
halos, in addition to the complications arising with the 
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high and low resolution simulations of set A, prompted 
us to conduct a second set of simulations. Set B con- 
sisted of 408 3 particles giving a coarse particle mass of 
3.8 x 10 8 M Q and a force resolution of 4.4 kpc. We ran 
HOP and identified halos with masses 0.8- 2.2 x 10 12 M o . 
For each halo we also found the nearest neighboring halo 
with mass > 0.8 x 1O 12 M . We selected a halo whose 
nearest massive neighbor was at least 5 Mpc away and 
visually verified the halo was indeed isolated. A rect- 
angular refinement level 6.1 x 7.0 x 7.9 Mpc was placed 
over this halo's particles in the initial conditions. Low 
and high resolutions were conducted with the same mass 
and force resolutions as set A. The low resolution sim- 
ulations used 16,515,072 (~ 255 3 ) particles in the re- 
finement level while high resolution used 132, 120, 576 
(~ 510 3 ) particles in the refinement level. 

A third set of CDM only, low resolution simulations 
was run to further explore the scatter between the sub- 
halo distributions of different realizations of Milky Way- 
type halos and to explore the possibility of a bias intro- 
duced by the use of the BBKS transfer function. Set C 
consisted of 408 3 particles but the CDM transfer function 
was generated from the LINGER software in GRAFIC2 
[30j after correcting a bug where the power spectrum for 
baryons was used for dark matter when calculating the 
transfer function. We used the AMIGA'S Halo Finder 
(AHF) software [66[ to find MW-sized halos with no equal 
sized neighbor within two virial radii (defined below). 
Nine halos were selected for refinement at low resolution 
from a variety of environments, low density with few large 
halos to high density with many large halos. The rectan- 
gular refinement regions had lengths 7.5 — 15.8 Mpc and 
31,752,192- 69,009,408 (316 3 - 410 3 ) particles. 

Table U summarizes the properties calculated by AHF 
for all of our simulated Milky Way halos at z = 0. We 
write i?A to mean the radius enclosing an overdensity A 
times the critical value. The mass and number of parti- 
cles inside Ra are Ma and Na, respectively; va is the 
circular velocity v A = GMa / Ra at Ra , and v max is the 
maximum circular velocity of the halo. We use the value 
A = 178f^ n 4 = 100 [13 (which is also very close to the 
value using the definition from [68| ) for the virial radius of 
our MWs and consider subhalos within i?ioo when com- 
paring to other published work. The mass, radius, and 
velocity at A = 50 are also used in the literature and 
these values are also listed in Table HI 

Figure [3] shows the density profiles of the A and B 
Milky Way halos calculated by breaking the halos into 
spherical shells. Small differences between the high and 
low resolution set A halos caused by the merging neighbor 
are apparent but generally the profiles are very similar 
across all simulations of each set. We do not see an inner 
flattening of the halos in the WDM simulations because 
we did not add thermal motions to our particles. Fig- 
ure 2] shows portraits of the Milky Way in our set B high 
resolution simulations. 
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FIG. 3. Density profile of Milky Way halos in the set A and 
set B CDM and WDM simulations. Thick lines are the high 
resolution simulations. Set B simulations are at top, the set 
A and the WDM cosmologies in each set have been offset 
downward for clarity. The profiles areplotted starting from 
the convergence radius of Power et al. [6!| for both resolutions 
(vertical lines). 



A. Identification of Satellites 

We used the AHF halo finding software [66| to find the 
gravitationally bound dark matter halos in our A-body 
simulations. Unbound particles were iteratively removed 
and we selected gravitationally bound halos with ten or 
more particles. 

AHF calculates properties of the halos it finds such as 
the total mass and the maximum circular velocity. For 
our purposes the maximum circular velocity is a better 
characteristic of a subhalo than the mass because quan- 
tifying the outer boundary of a subhalo embedded in a 
larger halo is somewhat arbitrary and can introduce sys- 
tematic errors. The maximum circular velocity however 
typically occurs at a radius well inside the subhalo out- 
skirts. 



III. RESULTS 

A. Satellite Distribution Functions 

We first compare our CDM simulations to other CDM 
simulations in the literature. Figure [S] shows the cumula- 
tive mass functions, N{> M su b), for subhalos within R§q 
for the set A and B MWs. Poisson statistic error bars 
have been added to the high resolution simulations and 
fit by N cx M~P. The values of 6 (0.9 and 0.95) agree 
with other published work that find values of 0.7 — 1.0 
0, 170l - l77j . At both high and low resolution the simulated 
mass functions turn away from the fit at masses below 
about 200 times the mass resolution of the simulation; 
this was also seen in the Via Lactea simulation [75[ • Also 



5 



TABLE I. Properties of simulated Milky Way halos. 



Simulation 
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1 keV lo 


1.4983 
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1.6615 


377.18 


137.64 


180.04 


2 
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2 


264 
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309.49 
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22 
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4 keV hi 


1.8261 
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2.0383 


403.77 


147.34 


189.69 
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22 
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19 
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22 
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1 keV hi 
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309.33 


2.0244 
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20 
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22 
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CDM lo 


1.9005 


312.84 
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2 
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2 
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plotted are the mass functions for the set C simulations. 
The subhalo abundances of the set A and B halos are 
within the halo to halo scatter and are consistent with 
the set C simulations. 

The cumulative maximum circular velocity functions 
for subhalos within i?ioo are plotted in Figure [5] The 
maximum velocities of the subhalos have been normal- 
ized by the maximum circular velocity of their host MW. 
The shaded region shows the minimum and maximum 
(lighter) and ±1ct (darker) from the mean of the 68 ha- 
los with masses 1.5 — 3 x 1Q 12 Mq in the simulation of 
Ishiyama et al. (78|. We use the fit to the density pro- 
file of the Via Lactea II halo [29[ to estimate its i?ioo 
(298 kpc) and from the published subhalo catalog we cal- 



culate and plot the Via Lactea II velocity function as the 
dashed line. The solid straight line is the fitting formula 
from the Bolshoi simulation [79l |. 

iV(> z) = 1.7 x 10-y^ osr r^ (9) 

= Vjnax I Vrnax.hosti (^^) 

applied to our high resolution halos which provides an 
excellent fit (the difference between the fit for the set A 
and B v max ,host is less than the thickness of the line). 
Via Lactea II used the same cosmological parameters as 
our simulations and their subhalo abundance agrees well 
with our simulations. Our simulations arc consistent with 
the Ishiyama et al. simulation but are systematically on 
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FIG. 4. Portraits of the simulated Milky Way halo at z = in 
of the MW centers is shown. 

the low end of their distribution. This is likely due to the 
different cosmology used in the Ishiyama et al. simulation 
(discussed below). 

Figure [7] also plots the cumulative velocity functions 
but includes all subhalos within i? 50 and the subhalo ve- 
locities have been normalized by the circular velocity at 
i?50 of their host MW. The Ishiyama et al. halos are 
again plotted as in Figure [5] as well as Via Lactea II. 
The solid straight line is the result from the Aquarius 
simulations [73] ■ Again there is good agreement between 
our simulations and Via Lactea II but an offset between 
our simulations and those of Ishiyama et al. and Aquar- 
ius. To first order, the abundance of halos of any size 
depends on the power spectrum of density perturbations 
which depends on the normalization, erg, and the tilt of 
the power spectrum, n s . Larger values of either param- 
eter increases the power on small scales and leads to a 
larger number of satellies for a given mass and v max of 
the host. The values (as = 0.9, n s = 1) were used in the 
Aquarius simulations and (0.8, 1) were used by Ishiyama 




the set B high resolution simulations. Structure within 500 kpc 



et al. Both are significantly greater than our adopted 
values of (0.74, 0.951), and this is the likely cause of the 
abundance offset. 

We adopted a WMAP3 cosmology to facilitate com- 
parison to the Via Lactea II simulation. The WMAP3 
values of n s , as, and fi m are 1.0, 2.9, and 2.5 standard 
deviations below the latest WMAP7 values HfJ. The 
Bolshoi simulation used parameters in agreement with 
WMAP7 and constraints from other cosmology projects. 
A comparison of the subhalo abundances of 4960 Bolshoi 
halos with circular velocities and masses comparable to 
the Via Lactea II halo indicated Bolshoi has more sub- 
halos by about 10%. Although Via Lactea II is just one 
halo and may not be representative of the average for a 
WMAP3 cosmology, this agrees with expectations from 
the 10% smaller value of as used by Via Lactea II. We 
used the same value of as as Via Lactea II but the Bolshoi 
fitting formula applied to our high resolution simulations 
in Figure [5] provides an excellent fit with no indication 
of an offset. This could be because, as we show in the 
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FIG. 5. Cumulative mass functions for subhalos within i?so for the CDM set A (left) and set B (right) simulations. Subhalo 
masses have been normalized by the M50 mass of the host. Poisson error bars have been added to the high resolution simulations 
and fit with a straight line. Both high and low resolution (dotted lines) turn away from the straight fit below about 200 times 
the mass resolutions of the simulations (short vertical lines). Mass functions for the set C halos have also been added (thin 
lines) and show the A and B abundances are within the halo to halo scatter. 
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FIG. 6. Cumulative velocity functions for subhalos within 
i?ioo in our low resolution set C (thin lines) and high resolu- 
tion set A and set B (thick lines) MW halos. Subhalo circular 
velocities have been normalized to the maximum circular ve- 
locity of the host halo. The dashed line is the subhalo velocity 
function of Via Lactea II, the straight solid line is the fitting 
formula from the Bolshoi simulation applied to our A and B 
halos. The shaded regions show the minimum and maximum 
and ±lcr from the mean of the 68 MW-sized halos of Ishiyama 
et al. [z3 
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FIG. 7. Cumulative velocity functions for subhalos within 
Rso in our low resolution set C (thin lines) and high resolu- 
tion set A and set B (thick lines) simulations. Subhalo circu- 
lar velocities have been normalized to the circular velocity of 
the host halo at a radius enclosing an overdensity of A = 50. 
The dashed line is the velocity function for Via Lactea II 
and the straight solid line is the average abundance from the 
Aquarius simulations [7^| . The shaded region is the min- 
imum/maximum range and ±I<r about the mean for halos 
from the simulations of Ishiyama et al. [7S| . 



Appendix, the BBKS power spectrum used in our high 
resolution simulations has about 10% more power at sub- 
Galactic scales. Below we argue that an intrinsic scatter 
in subhalo abundance of 30% (la) is reasonable to adopt 
and we conclude this can account for variations in the 
adopted cosmology without the need for a separate cor- 
rection. 



From our simulations and those of Ishiyama et al. in 
Figures [5] and [7] it is clear that the subhalo abundances of 
halos of similar sizes have a scatter. For a given cosmol- 
ogy the scatter in abundance includes an intrinsic scat- 
ter and a statistical scatter from the number of subhalos. 
The Aquarius simulation suite included 6 MW-sized ha- 
los simulated at very high resolution. For subhalos within 
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i?50j a t high values of the abundances where the statis- 
tical scatter is small, the la intrinsic scatter was deter- 
mined to be 10%. In Figurc[7]the scatter in the Ishiyama 
et al. abundances decreases for increasing N and appears 
to be converging to the 10% found in Aquarius. However 
the variation in the Ishiyama et al. abundances in Fig- 
ured)] is clearly converging to a larger value. As argued 
in [78| , the smaller scatter in Figure [7] can be explained 
by the inclusion of subhalos at distances up to R$o which 
are outside the virial radius and, hence, their evolution 
has not been affected by the structure of the host halo. 
Using V50 to normalize the subhalo velocities can also re- 
duce scatter since, unlike v ma x, it is less dependent on 
the central concentration. We will be interested in the 
number of subhalos in the inner regions of the host MW 
which are expected to be sensitive to the host concentra- 
tion so we will adopt the higher value for the intrinsic 
variation in the number of subhalo from Figure [5] which 
we estimate to be about 30% (la) after subtracting the 
Poissonian statistical scatter expected from the number 
of subhalos. 

In Figure [8] the cumulative circular velocity functions 
for subhalos within i?ioo for the high resolution set A 
and set B CDM and WDM simulations are plotted. The 
set A abundances have been increased 7% to normalize 
the CDM abundances to those of the set B simulation 
and illustrate that the relative suppression of subhalo 
abundances for each WDM simulation compared to CDM 
is the same across both simulation sets. The straight 
line is the Bolshoi fitting function applied to the set B 
CDM halo. The vertical lines in Figure [5] show where 
Vmax = 6 and 8 km/s. Below 8 km/s the high resolu- 
tion CDM simulations begin to fall away from the Bol- 
shoi line due to the resolution limits of the simulations. 
For v max > 8 km/s our simulations are reasonably com- 
plete within i?ioo of each Milky Way although numerical 
destruction of a small fraction of satellites in the inner 
Milky Way would not be apparent in Figure especially 
for the CDM and 4 keV cosmologies. Before comparing 
our simulations to observations we need to determine to 
what distance our simulations are convergent. 



B. Convergence Study 

Satellites orbiting in the halo of a larger galaxy are de- 
stroyed by tidal stripping and heating through encounters 
with other satellites. Satellites in simulations are also de- 
stroyed artificially by numerical effects that become dom- 
inant for poorly resolved halos in the inner halo region. 
There will therefore be a radius inside of which our sim- 
ulations will not converge to a realistic representation of 
the actual Milky Way. 

To determine the convergence of our simulations and 
have an idea of the variance of the results, we performed 
simulations at lower and higher resolution of two differ- 
ent realizations of a Milky Way-sized galaxy. We per- 
formed convergence studies following the argument elu- 
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FIG. 8. Cumulative velocity functions for satellites in our high 
resolution set A and set B CDM and WDM simulations. The 
set A abundances (thin lines) have been increased by 7% to 
normalize the CDM abundances to those of the set B and show 
that the relative suppression of halos in WDM cosmologies 
compared to CDM is similar across both simulations. The 
straight gray line is the Bolshoi fitting formula applied to the 
set B halo. 

cidatcd below, in combination with results of published 
high resolution simulations found in the literature. Using 
the work of Moore et al., De Lucia et al., and Ishiyama et 
al. 0, [zl, [78| , we will assume that the shape of the cumu- 
lative satellite velocity function for host halos of different 
masses is nearly constant and the total number of satel- 
lites scales linearly with the host mass. If the simulations 
are convergent, the cumulative circular velocity function 
for satellites, N(R), within a given Galactocentric radius, 
R, should be proportional to the enclosed mass, M(R), 
and a function of R that represents the fraction of satel- 
lites that survive destruction from physical effects: 

N(R) = f(R)M(R), (ff) 

where f(R) oc R a . The normalization of f(R) can be set 
using values of N(R) and M{R) at a distance Ro: 

The velocity functions normalized in this way will be con- 
stant with radius where the simulations are convergent. 
Where numerical effects destroy satellites the velocity 
functions will normalize to a lower value. We expect 
that a is constant because there is no characteristic scale 
for the destruction rate in dark matter only simulations. 
Hence, a can be determined at large radii where conver- 
gence is certain. 

Figure [9] shows the normalized velocity functions for 
the simulations. The normalization constants M ,R 
have been chosen at 200 kpc and the value of a (0.55) 
was adjusted by hand until a good fit was achieved for 
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the set B velocity functions above 200 kpc in the high res- 
olution CDM cosmology at circular velocities > 6 km/s 
(vertical lines). This a also provides a good fit for the 
WDM cosmologies and for the set A simulations, al- 
though the 1 and 2 keV velocity functions have a wider 
scatter due to the smaller numbers of satellites in these 
simulations. The tuwdm = 4 keV simulation is conver- 
gent for v max > 6 km/s to distances > 100 kpc. At 
75 kpc the effects of numerical resolution are apparent. 
The same value of a has been used in the normaliza- 
tion of the low resolution sets and appears to provide a 
good fit for the velocity functions > 200 — 250 kpc. The 
effects of numerical resolution on the destruction of satel- 
lites are apparent at larger distances in these simulations: 
< 200 kpc for CDM and < 150 kpc for WDM. 

The 1 and 2 keV WDM velocity functions in Figures [8] 
and [5] show a flattening when going from high to low 
velocities until about 6 km/s, below which the number 
of subhalos increases gre atly . This is a generic feature 
of WDM simulations [l(J, HH and is usually explained as 
top-down fragmentation of matter filaments. Given that 
WDM simulations are known to form nonphysical halos 
along filaments [24|,[25[ , it is likely the low velocity upturn 
in the velocity function of subhalos is actually caused by 
these numerical artifacts accreting onto the MW. Since 
these nonphysical halos form with separations typical of 
the mean particle distances in the initial conditions, the 
number of halos should increase with the mass resolution 
of the simulation. Our low resolution simulations do not 
show clear evidence of upturns in the velocity functions, 
we require simulations with resolution higher than our 
high resolution sets to confirm this effect. Regardless, 
we consider only satellites with velocities greater than 
6 km/s in our high resolution simulations when deriving 
our constraints on the dark matter particle mass. 



C. Comparison to Observations 

Before the Sloan Digital Sky Survey there were only 
12 classically known satellite galaxies to the Milky Way. 
Sixteen new satellites have been discovered in the SDSS, 
currently in Data Release 7. We list all known Milky 
Way satellites in Table [TTJ We use the satellite distances 
given in Table HT1 as their Galactocentric distances. When 
comparing the observed satellites to the simulations, we 
must correct the SDSS dwarfs for completeness. The 
primary incompleteness of the SDSS is its sky coverage, 
28.3% (11663 deg 2 ). Second, being a magnitude limited 
survey, the SDSS has a luminosity bias. The detection ef- 
ficiency of dwarfs in the SDSS is a function of dwarf size, 
luminosity, distance, and Galactic latitude as shown by 
Walsh et al. 18 2| . An approximate expression is given in 
Tollerud et al. j83| (using the work of Koposov et al. |84j |) 
for the distance which galaxies of luminosity > L are 
completely detected: d « 66kpc(L/1000 Lq) 1 ^ 2 . Galax- 
ies with L > 10 4 Lq should be approximately complete 
to 200 kpc, with L > 2300 L to 100 kpc. The distance 



range 100—200 kpc is thus suited for comparisons because 
our simulations are convergent and the observations are 
nearly, but not quite, complete. For our analysis we will 
only use satellites with distances < 200 kpc. For compar- 



TABLE II. Summary of known Milky Way satellites. 



Name 


dist 
[kpc] 


O star 

[km/s] 


My 


References 


Classical (pre-SDSS) 


Sagittar 


24 ±2 


11.4 ±0.7 


-13.4 


[85] 


LMC 


49 ± 2 




-18.4 


[85] 


SMC 


58 ±2 




-17.0 


[85] 


Ursa Minor 


66 ±3 


9.3 ± 1.8 


-8.9 


[85] 


Draco 


79 ±4 


9.5 ± 1.6 


-8.8 


[85] 


Sculptor 


79 ±4 


6.6 ±0.7 


-11.1 


[85] 


Sextans 


86 ±4 


6.6 ±0.7 


-9.5 


[85] 


Carina 


94 ±5 


6.8 ± 1.6 


-9.3 


[85] 


Fornax 


138 ±8 


10.5 ±1.5 


-13.2 


[85] 


Leo II 


205 ± 12 


6.7 ±1.1 


-9.6 


[85] 


Leo I 


270 ± 30 


8.8 ±0.9 


-11.9 


[85] 


Phoenix 


405 ± 15 




-10.1 


[85] 


SDSS discovered 


Segue I 


23 ±2 


4.3 ± 1.2 


-1.5 


[86] 


Ursa Major II 


30 ±5 


6.7 ± 1.4 


-3.8 


[87, 88] 


Segue II 


~ 35 


3.4 ±2.0 


-2.5 


[89J 


Willman I 


38 ±7 




-2.5 


[87, 90] 


Coma Berenics 


44 ±4 


4.6 ±0.8 


-3.7 


[88, 9^ 


Bootes II 


60 ± 10 




-3.1 


[92] 


Bootes I 


62 ±3 


6.5l?;° 


-5.8 


[87J 


Pisces I 


80 ± 14 






[93, 94] 


Ursa Major I 


10611 


7.6 ± 1.0 


-5.6 


[88] 


Hercules 


140±il 


5.1 ±0.9 


-6.0 


[88, 9^ 


Canes Venatici II 1501^ 


4.6 ± 1.0 


-4.8 


[88, 91] 


Leo IV 


i6o±it 


3.3 ± 1.7 


-5.8 


[88, 91] 


Leo V 


175 ±9 


2.4 ± 1.8 


-5.2 


[95, 96] 


Pisces II 


~ 180 




-5.0 


[97] 


Canes Venatici I 


220t^ 


7.6 ±0.4 


-7.9 


[88, 98] 


Leo T 


~ 420 


7.5 ± 1.6 


-7.1 


[88, 99] 



ison with the simulations, we do a first order correction 
for the sky coverage of the SDSS assuming an isotropic 
distribution of satellites. The SDSS covers 0.283 of the 
sky, so for every satellite detected we assume there are 
actually 3.54 satellites at that distance. Combined with 
the classic Milky Way satellites this forms our observed 
data set. We also implemented a conservative luminos- 
ity correction for the SDSS dwarfs using the formulas 
in Walsh et al. [HI]. This adds two satellites making a 
total of 61 ± 13 satellites within 200 kpc where the er- 
ror only includes Poisson statistics of the SDSS dwarfs: 
3.5Ay/Ngjjgg. The extra two satellites are at distances 
150 — 200 kpc and do not affect our conclusions. We note 



10 



Set B Hi-res CDM 

5keV 





Set A Hi-res CDM 

5keV 
4keV 



5 10 20 

»™, [km/s] 




"ma* [km/s] 



FIG. 9. (left) Velocity functions for set B high (top) and low {bottom) resolution simulations normalized with Eq. (|12[) . Solid 
lines are R = 400, 300, 250, and 200 kpc (thick), dotted line is R = 150 kpc, dashed line is R = 100 kpc, dot-dashed line is 
R = 75 kpc. The WDM cosmologies have been shifted down vertically for clarity. The value a = 0.55 was set by the high 
resolution simulation and provides good normalization for the low resolution as well but the effects of incompleteness become 
apparent at much larger radii (150 — 200 kpc compared to 75 — 100 kpc for high resolution), (right) Same as the left panel, but 
for the set A high (top) and low (bottom) resolution simulations. The value a = 0.55 also provides good normalization for this 
set of halos. 



that we get the same result if we only include dwarfs 
detected in the Data Release 5 footprint and correct for 
the smaller sky coverage (8000 dcg 2 ). Wc also note the 
formulas in Walsh et al. [82[ assume the size-luminosity 
distribution of known dwarfs is representative of all satel- 
lites. There may be a population of dwarfs with s urfac e 
brig htnesses below the detection limit of the SDSS flOOr - 

Willman 1 is an exceptional case in that it may not 
be a dark matter dominated dwarf galaxy but a globular 
cluster undergoing tidal disruption. Its stellar velocity 
dispersion implies a large mass to light ratio like other 
dwarf sphcroidals and it has a size and luminosity inter- 
mediate between MW dwarfs and globular clusters [90j |. 
but unresolved binaries and tidal heating may contami- 
nate the velocity dispersion and lead to an overestimated 
mass. Although it has a large mctallicity spread unlike 
the stellar popul ation of a globular cluster [83], follow- 
up spectroscopy [104| ] suggests there may be contamina- 
tion by foreground stars and when these are excluded 
the metallicity spread can be consistent with a metal- 
poor globular cluster. The detection of an X-ray emission 
line fro m decayin g dark matter in Willm an 1 has been 
claimed [l05l 1 106| but is provisional [l07[ | . If confirmed 



this would indicate a dark matter halo in Willman 1 and 
confirm its status as a MW satellite. When deriving con- 
straints on the dark matter particle mass we will consider 
both including and excluding Willman 1 as a Milky Way 
satellite. 

When comparing observations and simulations we will 
apply cuts to the simulated subhalos and consider only 
those with velocites above 6 and 8 km/s. As discussed 
in the previous section this is to avoid potential contam- 
ination from numerical effects in the WDM simulations. 
That these velocity cuts are a reasonable estimate of the 
minimum v max of the dark matter halos the observed 
galaxies are presumably e mbed ded in can be shown as fol- 
lows. Ricotti and Gnedin [lOOj found in simulations that 
the maximum circular velocities of satellites are at least 
twice the velocity dispersion of the stellar component. 
Assuming the stellar velocity dispersions of the observed 
dwarfs are times the line-of-sight velocity dispersions 
(cstor m Table [TTJ) , then all dwarfs with measured veloc- 
ity dispersions have maximum circular velocities greater 
than 8 km/s. We assume dwarfs without measured veloc- 
ity dispersions are similar to the other known dwarfs and 
conservatively conclude all dwarfs reside in dark matter 
halos with v max greater than 8 km/s. An alternative ap- 
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proach is the work of Wolf et al. |l08| relating the circular 
velocity at half light radius to (J star' v c (r 1 / 2 ) = V3a star - 
All the observed dwarfs except Leo V have circular veloc- 
ities at half light radius about 6 km/s or greater. Since 
the maximum circular velocity must be greater than or 
equal to the half light circular velocity, it is also reason- 
able to consider that all observed dwarfs reside in halos 
with v max > 6 km/s. We stress that these cuts reflect 
the need to reduce the numerical effects of the nonphys- 
ical halos in WDM simulations that dominate our high 
resolution simulations at subhalo velocities below 6 km/s 
rather than an assumption on the relationship between 
luminous satellites and dark matter halos. 

In the left panel of Figure [TU] we plot histograms of 
the number of satellites with distance for the observed 
and simulation data sets with a 6 km/s maximum cir- 
cular velocity cut. The upward arrows on the observed 
data bars indicate these are only lower limits due to the 
surface brightness limits of the SDSS; it is possible there 
are more dwarfs yet to be discovered. Willman 1 has 
not been included as a MW satellite in this hgurc. The 
6 km/s cut to the simulation data assures the high reso- 
lution simulations are convergent to at least r = 100 kpc. 
The low resolution simulations are also plotted in these 
figures but they are convergent only to 150 kpc. Focusing 
on the 100 — 200 kpc bins it is clear the 1 keV has far 
too few satellites to match the observations. The 2 keV 
simulations are consistent with the observations in the 
100 — 150 kpc bin but to be generally consistent with the 
observations would require the simulations to be signifi- 
cantly incomplete below 100 kpc or the sky correction to 
have overestimated the number of satellites in the inner 
Milky Way. The 4 keV simulations can be consistent with 
the observations although they may require some of the 
dark matter halos to not host luminous galaxies. Strong 
conclusions cannot be drawn from this plot because it is 
not clear how variance in the satellite abundances for the 
simulated halos may affect the results. 

The number of satellites in the simulations can be cor- 
rected for completeness using the convergence equation, 
Eq. (fT"2"|) . We used the mass and number of satellites 
inside R§q for the normalization and calculated the num- 
ber of satellites in 50 kpc bins for the high resolution 
simulations. The results are shown in the right panel of 
Figure QJJ The results are very similar across cosmolo- 
gies, a nearly constant number of satellites per bin from 
50 — 200 kpc with about half as many in the — 50 kpc 
bin. The plots are also very similar to the simulation 
data plots except for the — 50 kpc bin in the 4 keV 
cosmology where about ten satellites were destroyed by 
numerical effects in the set B halo. The — 50 kpc bin is 
most important for constraining the dark matter particle 
mass because the observations are most complete in this 
bin. 

In the top panel of Figure [TTJ we plot the number of 
satellites in the — 50 kpc bin calculated from the conver- 
gence equation for the high resolution set B Milky Way 
as a function of the particle mass with a 6 km/s veloc- 



ity cut to the subhalos. We set the variance in satellite 
abundances equal to that for a 30% intrinsic rms scat- 
ter plus a Poissonian variance in the number of satel- 
lites. The dark and light shaded regions in the plot show 
the la and 2a ranges, respectively. The solid horizontal 
line shows the number of observed satellites after apply- 
ing sky coverage corrections to the SDSS data (excluding 
Willman 1) while the dashed and dotted horizontal lines 
show the 68% and 95% confidence lower limits from Pois- 
son statistics. In the bottom panel of Figure [TT] we plot 
the difference in the number of observed and simulated 
satellites with la and 2a limits (shaded regions) from 
combining the variances of the observed and simulated 
data sets. The number of satellites in simulation must 
be at least equal to the number of observed satellites, 
therefore where this quantity equals zero defines a lower 
limit on the dark matter particle mass. The arrowed lines 
indicate the lower limits at la and 2a for this case of the 
set B simulation with a 6 km/s cut to the subhalos and 
excluding Willman 1 from the observed set. 

We repeated the same analysis using a v max > 8 km/s 
cut to the simulation data. We also considered the effects 
when Willman 1 was included in the observed data set 
for both the 6 km/s and 8 km/s analysis. We present 
the results in Tabic Mil The set B halo is slightly 
more abundant in subhalos but both the set A and B 
simulations give the same results to within about 10%. 
Rather than take the average of the two simulations we 
will simply adopt the more conservative of the two con- 
straints. In the most conservative case, where Willman 1 

TABLE III. Dark matter particle mass constraints (in keV) 
from the high resolution set A and B MW halos. Constraints 
for simulated subhalo v max cuts of 8 and 6 km/s and including 
or excluding Willman 1 from the observed data set are given. 





v m ax > 8 km/s 


Umax ^ 


6 km/s 


Will 1? 


Included 


Excluded 


Included 


Excluded 


MW 


A B 


A B 


A B 


A B 


2a 


> 3.6 > 3.3 


> 2.9 > 2.7 


> 3.0 > 2.6 


> 2.5 > 2.3 


la 


> 5 > 4.6 


> 4.2 > 3.8 


> 3.9 > 3.4 


> 3.3 > 2.9 



is not a dark matter dominated dwarf galaxy and all ob- 
served satellites correspond to dark matter halos with 
v-max > 6 km/s, we can say ttiwdm > 2.3 keV with 95% 
confidence. We adopt this as our formal limit for this 
work. 



IV. DISCUSSION 

We found that a model with m w d m = 4 keV produces 
the best fit to observations at < 50 kpc, i.e. this model 
has a number of dark matter satellites equal to the num- 
ber of observed luminous satellites. However, due to 
the large uncertainties in the number of observed satel- 
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FIG. 10. (left) Number of satellites with distance from the Milky Way grouped into 50 kpc bins. The simulation satellites 
have been cut by circular velocity > 6 km/s. Solid lines show the data from the high resolution and dashed lines show the low 
resolution simulations, thick lines are the set B simulations. The bars with arrows are the observed satellites after correcting 
for the sky coverage of the SDSS but not including Willman 1. Observations are incomplete at distances greater than about 
50 — 100 kpc (depending on the luminosity and surface brightness of the dwarf), while simulations have not converged for less 
than about 100 kpc for high resolution and 150 kpc for low resolution, (right) Number of satellites with distance from the Milky 
Way like the plot at left but calculated from the convergence equation [Eq. [T2] . The convergence correction affects mainly the 
4 keV 0-50 kpc bin. 



lites due to partial sky coverage and on the number of 
simulated satellites due to Poisson and intrinsic scatter, 
that partially reflects observational uncertainties on the 
mass and v max of the Milky Way, we find much weaker 
lower limits on m w d m than 4 keV. In the future however, 
the lower limit on m w d m will improve as observations of 
MW satellites become more complete. The scatter of the 
simulation can also be reduced using constrained simu- 
lations of the Local Group (also including the effect of 
baryons) in combination with more accurate determina- 
tion of the mass, rotation curve, and concentration of the 
Milky Way. 

Considering the various uncertainties in the number of 
observed and simulated satellites, we found a conserva- 
tive lower limit of ttiwdm > 2.3 keV (2a) on the dark 
matter particle mass. We also found the 1 keV WDM 
simulations have too few satellites to match the Milky 
Way observations. This agrees with the semianalytic 
modeling and Milky Way satellite luminosity functions 
in WDM cosmologies work of Maccio and Fontanot 14| ; 
however, we only apply a cut to the simulated halos to 
avoid numerical effects and do not make assumptions on 
how the dark matter halos arc populated by luminous 
galaxies. 

Our result can also be compared to limits on the 
particle mass from the Lyman-a forest in high redshift 
quasars. Lyman-a absorption by neutral hydrogen along 
the line of sight to distant quasars over redshifts 2-6 
probes the matter power spectrum in the mildl y nonlin- 
car regime on scales 1-80 Mpc/h. Vicl et al. f5^. fl09l [TTo| 
have numerically modeled the Lyman-a forest flux power 



spectra for varied cosmological parameters and compared 
to observed quasar forests to obtain lower lim its o n the 
dark matter particle mass. Their 2006 work jl09j used 
low resolution s pect ra for 3035 quasars (2.2 < z < 4.2) 
from the SDSS 11 if and found a 2a lower limit of 2 keV 
for a thermal WDM particle. This limit agrees with our 
results that a 2.3 keV particle is the lower limit that can 
reproduce the observed number of Milky Way sate llites 
and agrees with the Lyman-a work of Scljak et al. |l!2| 
who find a 2a limit > 2.5 ke V fo r a thermal particle. 
The latest work of Viel et al. jllOj uses high resolution 
spectra for 55 quasars (2.0 < z < 6.4) from the Keck 
HIRES spectrograph in addition to the SDSS quasars. 
With the new data they report a l ower limit of 4 keV 
(2a). A caveat arises in Viel et al. |l!3| . who show the 
flux power spectrum from the SDSS data prefer larger 
values of the intcrgalactic medium (IGM) temperature at 
mean density than expected from photoionization. The 
flux power spectrum temperature is also higher than that 
derived from an analysis of the flux probability distri- 
bution function of 18 high resolution spectra from the 
Very Large Telescope and also higher than constraints 
from the widths of thermally broadened absorption lines 
Il4lll5| . This could be explained by an unaccounted for 



systematic error in the SDSS flux power spectrum data 
which may also affect the derived dark matter particle 
mass limits. 

Using the scaling relation for sterile neutrinos we find a 
lower limit m s > 13.3 keV with 95% confidence for a DW 
produced sterile neutrino particle. Scaling to the other 
production mechanisms we get m s > 8.9 keV for the SF 
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Dark matter particle mass [keV] 

FIG. 11. top Number of satellites from 0-50 kpc calculated 
from the convergence equation in the set B high resolution 
simulation for v max > 6 km/s (sloped solid line) with la and 
2a limits (shaded regions) compared to the observed number 
of Milky Way satellites, excluding Willman 1, with correction 
for sky coverage (solid horizontal line) and at 68% and 95% 
confidence (dashed and dotted horizontal lines), bottom The 
observed number of satellites minus the number of satellites 
in simulation with la and 2a limits (shaded regions). The 
number of satellites in simulation must be greater than or 
equal to the number of observed satellites, therefore where 
this quantity equals zero sets a lower limit to the dark matter 
particle mass (arrows). 



ergy E 1 — m s /2. X-ray observations of the diffuse X- 

ray background Ill6l and dark matter halos in clusters 
[117141201 . M31 [ml, dwarf spheroidal galaxies [l22r - 
Il25| . and the halo of the Milky Way HH [lH, [l27| have 
all been used to set constraints on the sterile neutrino 
mass. Observatio ns of the diffuse X-ray background have 
set m s < 9.3 keV [l!6j . while the Virgo and Com a clus- 
ters have been used to set m s < 6.3 keV [l 1 7l ] which 
also agrees with limi ts fro m the bullet cluster, IE 0657- 
56, m s < 6.3 keV [l20t ] and is cl ose t o results from 
the Milky Way halo m s < 5.7 keV fl27j . Tighter con- 
straints have been determined from M31 observations 
m s < 3.5 keV |l2l| a nd f rom the dwarf spheroidal Ursa 
Minor m s < 2.5 keV [l25[. 

These upper limits are well below the lower limits de- 
rived in this work and from Lyman-a observations and 
seem to rule out the DW and SF production mecha- 
nisms. However, all of these mass limits, including the 
constraints set in this work, arc model dependent and 
make certain assumptions. In general X-ray constraints 
depend on the sterile neutrino mass, the mixing angle 
with active neutrinos 9, and the cosmic matter density 
of sterile neutrinos f2 s . There are also assumptions about 
the initial conditions, that there were no sterile neutri- 
nos in the early Universe at temperatures > 1 GeV, there 
was no entropy dilution after creation, and no coupling 
to other particles. There are also uncertainties with the 
calculation of production rates because these occur at 
temperatures where the plasm a is neithe r well described 
by hadronic nor quark models [l!6l Il28l | . Depending on 
the assumptions made and the adopted production model 
the relationship between m s , (9, and £l s changes so that 
robust constraints cannot be placed on any one model 
parameter. 

There has also been a report of a detection of a dark 
matter X-ray emission lin e from W illman 1 consistent 
with m, = 5.0 ± 0.2 keV HM [Tol. This detection is 



provisional 1071 ] but if confirmed the limits derived in 
this work imply the sterile neutrinos are not produced 
by the DW mechanism or do not constitute the entirety 
of the dark matter. 



mechanism and m s > 3.0 keV for Higgs decay sterile 
neutrinos; however, we note this is not based on transfer 
function calculations for the SF and Higgs mechanisms 
but assumes a simple scaling for the average momen- 
tum for the different production mechanisms [34| . The 
Lyman-a forest observations discussed above in the con- 
text of a thermal particle also set limits on t he s terile 
neutrino mass. The 2006 work of Viel et al. Il09l sets 
m s > 11 keV and is similar to the Seljak et al. |l!2ll limit 
m s > 14 keV. The 2008 work of Viel et al. [110| sets the 
highest limit of m s > 28 keV but is subject to the caveats 
mentioned above. 

Sterile neutrinos are expected to radiatively decay to 
a lighter mass neutrino and a X-ray photon with en- 



V. SUMMARY 

We conducted TV-body simulations of the formation of 
a MW-sized dark matter halo in CDM and WDM cos- 
mologies. Such simulations are complicated by the for- 
mation of nonphysical small mass halos due to the dis- 
creteness of the initial conditions but with sufficient res- 
olution they are only important at small scales and can 
be avoided with an appropriate circular velocity cut. 

We studied the number of satellite halos as a func- 
tion of distance from the MW. The 4 keV WDM sim- 
ulation can adequately reproduce the observed number 
of satellites at hundreds of kiloparsecs while the 2 keV 
simulation is slightly deficient and the 1 keV severely de- 
ficient. Our high resolution simulations followed the for- 
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mation of two MW-sized halos. Numerical simulations 
of MW-sized halos show significant variance in the num- 
ber of satellites, an effect that can be easily quantified 
using published studies and was incorporated in our re- 
sults. We calculated the number of satellites in the inner 
50 kpc, corrected for the effects of numerical destruction, 
and accounted for the variance by conservatively adopt- 
ing a 30% (ler) intrinsic scatter in the number of satellites 
in addition to a scatter from Poisson statistics. We also 
accounted for the uncertainty in the number of observed 
MW satellites due to the survey area of the SDSS and 
derived a very conservative lower limit on the dark mat- 
ter particle mass of > 2.3 keV (95% C.L.). This agrees 
wit h the earlier Lyman-a forest modeling work of Viel et 
al. [l09j | that mwDM > 2 keV but the two methods are 
independent and almost certainly are subject to differ- 
ent systematic errors if any exist. Their latest work [llOj ] 
raised the limit to tuwdm > 4 keV but problems with the 
derived IGM temperature and mean density may indicate 
the S PSS spectra they used may suffer a systematic error 

ma. 

Our lower limit of 2.3 keV for a thermal dark mat- 
ter particle scales to lower limits of 13.3, 8.9, 3.0 keV 
(95% C.L.) for DW, SF, and Higgs decay produced ster- 
ile neutrinos. Sterile neutrinos, if they exist, are expected 
to decay into X-rays and active neutrinos. Observations 
of the unresolved cosmic X-ray background and X-ray 
observations of dark matter halos on scales from dwarf 
galaxies to clusters set upper limits below our lower limit 
and the limits of Lyman-a forest modeling. These limits 
are derived under many assumptions and, in general, the 
constraints apply to a parameter space of m s , 9, and Q s . 

Our constraint is a conservative lower limit since we 
only correct the number of SDSS dwarfs for sky com- 
pleteness. An analysis that takes into account the sur- 
face brightness limits of the observational data may allow 
tighter constraints; however, the analysis would be some- 
what model dependent. We have also not included the ef- 



fects on subhalos of baryonic structures in the inner MW 
halo such as a disk. The presence of a disk could lead to 
greater subhalo destruction due to increased dynamical 
friction and tidal heating. By increasing the subhalo de- 
struction rate in the inner halo, disks would increase our 
lower bounds on the dark matter particle mass. The as- 
sumption of no disk is a conservative one and an analysis 
that includes a disk may allow tighter constraints. 

We have demonstrated how TV-body simulations of the 
MW and its satellites can set limits on the dark matter 
particle mass comparable to and independent of comple- 
mentary methods such as modeling the Lyman-a forest. 
These limits are helped greatly by the discovery of many 
new MW satellites in the SDSS. There may still be a pop- 
ulation of low luminosity, low surfac e br i ghtn ess dwarf 



galaxies undetectable by the SDSS |TooHlo3 [. Future 
surveys like LSST, DES, PanSTARRS, and SkyMapper 
have the potential to discover many more MW satellites 
and further improve constraints on the mass of the dark 
matter particle. Better constraints will result from the 
smaller uncertainty in the number of observed satellites 
achieved by improving the sky coverage and reducing lu- 
minosity corrections. In addition, the existence of a yet 
unknown population of even fainter satellites is not un- 
likely. 
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Appendix 

We used the BBKS formula for the CDM transfer func- 
tion when generating the initial conditions for our high 
resolution simulations. This formula assumes a baryon 
density of zero. Eisenstein and Hu [32j calculated trans- 
fer functions for CDM cosmologies that include baryon 
physics. We plot the power spectra for the fitting formula 
of Eisenstein and Hu and the spectrum calculated with 
the LINGER software (using Of, = 0.04) normalized by 
BBKS in Figurc[H With ft m = 0.238 a Milky Way-sized 
halo with mass ~ 2 x 10 12 M Q would form from a spher- 
ical region with diameter 4.8 Mpc (k = 0.28 h/Mpc); 
this is plotted along with the scale of the simulation box 
(90 Mpc) as solid vertical lines. Dashed vertical lines 
show the cell size in the refinement region of the low and 
high resolution simulations. Figure [T^] shows that, for a 
fixed value of ag, , BBKS underestimates power on scales 
k < 0.1 but the power spectra are nearly identical for 
scales < 14 Mpc with BBKS slightly power overabundant 
by ~ 10%. The set C halos showed subhalo abundance 
variations much greater than 10% and the BBKS power 
overabundance is much less than the 30% (lc) intrinsic 
scatter in subhalo abundance for MW-sized halos that 
we adopt. 

To check if BBKS might affect the number of satellites, 
we reran the set B high resolution CDM using initial con- 
ditions generated from the formula of Eisenstein and Hu. 
The panels of Figure [13] compare the velocity functions 
of satellites and show good agreement between the sim- 
ulations and within the scatter of the set C simulations. 
Based on this and the agreement between the BBKS and 
set C simulations seen in Sec. lIII Al we conclude that the 
use of BBKS has not introduced a systematic error into 
our results. 



19 




FIG. 12. Comparison of CDM power spectra calculated from the fitting formula of Eisenstein and Hu (EH97) and from the 
LINGER software normalized by BBKS. On scales k > 0.1 fe/Mpc (< 14 Mpc) the power spectra are nearly identical. The 
'MW vertical line is the diameter of a spherical region with density Q m pc enclosing a Milky Way-sized mass 2 x 10 12 Mq 
(~ 5 Mpc). This scale is well within the range where the power spectra are nearly equal. 




FIG. 13. Subhalo velocity function comparison for CDM high resolution set B simulations using fitting formula from BBKS 
and Eisenstein and Hu [thick lines) and the LINGER using set C simulations [thin lines), [left) Subhalos within -Rioo and 
velocities normalized by v ma x of their host. The straight sloped line is the fitting formula from the Bolshoi simulation, (right) 
Subhalos within Rso, normalized by vsq of their host. The subhalo abundances between the BBKS and EH97 simulations are 
in good agreement and within the scatter of the set C halos. 



